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ABSTRACT 

We investigate the formation and evolution of satellite galaxies using smoothed particle 
hydrodynamics (SPH) simulations of a Milky Way (MW)-like system, focussing on 
the best resolved examples, analogous to the classical MW satellites. Comparing with 
a pure dark matter simulation, we find that the condensation of baryons has had 
a relatively minor effect on the structure of the satellites' dark matter halos. The 
stellar mass that forms in each satellite agrees relatively well over three levels of 
resolution (a factor of ~ 64 in particle mass) and scales with (sub)halo mass in a 
similar way in an independent semi-analytical model. Our model provides a relatively 
good match to the average luminosity function of the MW and M31. To establish 
whether the potential wells of our satellites are realistic, we measure their masses 
within observationally determined half-light radii, finding that they have somewhat 
higher mass-to- light ratios than those derived for the MW dSphs from stellar kinematic 
data; the most massive examples are most discrepant. A statistical test yields a ~ 6 
percent probability that the simulated and observationally derived distributions of 
masses are consistent. If the satellite population of the MW is typical, our results 
could imply that feedback processes not properly captured by our simulations have 
reduced the central densities of subhalos, or that they initially formed with lower 
concentrations, as would be the case, for example, if the dark matter were made of 
warm, rather than cold particles. 
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1 INTRODUCTION 

Substantial progress has been made over the last few years 
in modelling the formation of galactic d ark matter halos us- 
ing high resolution N-body simulations |Springel et al.ll2008l : 
iDiemand etaH 120081 ; IStadel et al.1 bOOa T Hydrodynamical 
simulations of such systems inevitably lag behind in terms 
of resolution, but are now reaching a point where they can 
be used to investigate the detailed evolution of the baryonic 
component of satellite galaxies , as de monstrated by several 
recent studies. lOkamoto et al.l l|2010t) studied the effects of 
different feedback models on the chemical properties and lu- 
minosities of the satellit e populations arou n d thr ee Milky 
Way (MW)-mass halos. lOkamoto fc Frenkl (120091 ) showed 
that a combination of the early reionisation of pregalac- 
tic gas at high redshift and the injection of supernovae 
energy is sufficient to suppress star formation in the myr- 
iad of low mass subhalos that form in the ACDM cosmol- 
ogy, confirming results from earlier s e mi-analytic modelling 
(|Benson et alj2002l ; ISomerville!l2002D . IWadepuhl fc Springell 
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further argued that cosmic rays generated by super- 
novae may play an important role in suppressing star for- 
mation in satellites. 

From an observational point of view, t he release of data 
from the Sloan Digital Sky Survey (SDSS) (|York et alJl200Ch 
over the last decade has transformed the study of the Lo- 
cal Group satellites. The ~ 30 faint dwarf galaxies discov- 
ered using those data (e.g., Zuckcr ct al. 2004; Martin "et al.l 
l200fj ; iBelokurov et al.ll2007l ; iMcConnachie et al.ll200sh have 



prompted a new phase of detailed testing of current galaxy 
formation theories on smaller scales and in more detail than 
ever before. The SDSS data also reduced the discrepancy 
that existed between the number of low mass dark matter 
halos predicted by the ACDM cosmological model and the 
number of faint satellites identified around the MW: the 
'miss ing satellite problem ' l|Klvpfn et al.lll999l :l Moore et al.l 
1999). Over the same period, numerous theoretical mod - 
els fe.g.. lLi et aT1l2010l : lMacci6 et al.ll201ol : iFont et al.ll201ll ) 
have confirmed early conclusions that a combination of a 
photoionising background and feedback processes from su- 
pernovae (SNe) are capable of bringing the two into good 
agreement. 
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However, it is important to recognise that the satel- 
lite problem is not simply a statement that star formation 
must be suppressed in low mass halos. A more subtle, but 
equally important test of any cosmological model is whether 
the potential wells in which satellite galaxies form are capa- 
ble of suppor t ing st ellar systems with realistic kinematics. 
IStrigari et all (|201Ch demonstrated that all of the classical 
MW satellites for which high quality kinematic data are 
available are consistent with having formed in dark matter 
subhalos selected from the high resol ution ACDM N-bod y 
simulations of the Aquarius project (|Springel et al.| [2008). 
However, successful models must satisfy both constraints, 
producing realistic luminosity functions and forming stars 
in p otentials like those inf e rred f rom observations. 

iBovlan-Kolchin et alj (|20 1 if ) have recently argued that 
the most massive subhalos in high resolution simulations 
of cold dark matter halos are too concentrated to be able 
to host the bright est observed satellites of the Milky Way. 
lLovell et alj (|201 lh have shown that subhalos formed from 
warm, rather than cold, dark matter have suitably low con- 
centrations, but both they and Boylan-Kolchin et al. empha- 
sise that other solutions to the discrepancy are possible. One 
promising possibility is the mechanism originally proposed 
by Navarro, Eke & Frenk (1996a), whereby the rapid expul- 
sion of dense central gas in a starburst can unbind the inner 
parts of the halo, significantly reducing its concentration. 
We find an example of this process in one of the subhalos 
formed in the simulations analysed in this paper. 

From a theoretical perspective, ab initio hydrodynamic 
simulations are uniquely well-suited to investigating the ef- 
fects of galaxy formation on the dark matter halos of satel- 
lite galaxies. It has been known for some time that baryons 
may significantly alter the behaviour of dark matter on some 
scales. Dissipative processes such as gas cooling, star forma- 
tion and feedback decouple the dynamical evolution of the 
baryons from that of the dark matter. The resulting change 
in the shape of the overall potential in turn affects the phase 
space structure of the dark matter halo. 

Central concentrations of cold baryonic matter can in- 
duce an adiabatic, radial contraction of the cen tral re- 
gions (|Blumenthal et al.1 Il986l ; iGnedin et ail 120041 ). while 
the opposite effect can be achieved if dense clumps 
of baryonic material heat th e central distribut i on of 
dark matter (|Mo fc Mad 12004 iMashchenko et alj 120061 . 
2008), or if the blowout mechanism of Navarro, Eke and 
Frenk is effective (see also iGelato fc Sommer-Larsenlll999l; 



Table 1. Numerical parameters adopted for the three different 
resolution simulations: dark matter and gas particle masses and 
the gravitational softening in physical units. 



Gnedin fc Zhaol |2002l; 
20051 ; iGovernato et all 



Mo fc Mad 12004 iRead fc Gilmorel 
20101 ; iPontzen fc ; Governatd 1201^ 



This latter mechanism was originally proposed as a means 
to erase the central dark matter cusps in dwarf galax- 
ies, though whether or not such 'cored' profiles are re- 
quired by the observations remains a matter of ongo- 
ing debate IIGoerdt et al. 2006; Sanchez-Salcedo et al. 2006; 
IStrigari et al.ll2006l : lGilmore et al.ll2007l : I Walker et al.ll2009h . 
If baryons really do modify the dark matter in satellites 
on sub-kpc scales significantly, then the value of studying 
dwarf galaxies with post-processed N-body simulations may 
be very limited. 

In this paper, we make use of a model that has already 
had success in reproducing some properties of the Local 
Group satellites, including the shape and approximate nor- 
malisation of their luminosity function and the metallicity- 
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luminosity relation ( Ok amoto et al.l l2010t ) . In Section [2] we 
outline the details of our simulations, including the initial 
conditions, simulation code and modelling of various key 
baryonic physical processes. In Section [3] we expand on the 
theoretical predictions of the model by examining what ef- 
fect baryons have had on the dark matter profiles of satel- 
lites. We perform tests in Sectionf4]to ensure that key prop- 
erties of our satellite population do not depend on the nu- 
merical resolution. In Section [S] we compare several observ- 
able and derived properties of the simulated satellites to 
Local Group data. Finally, in Section [6] we discuss the evo- 
lution of one particularly interesting satellite in the simula- 
tion, which is dominated by its stellar component at z — 0. 
Our main results are summarised in Section [7] 



2 THE SIMULATIONS 

To investigate the properties of a simulated MW-satellite 
system, we select one of the six halo s from the Aquar- 
ius project described in Spr ingel et al.l (|2008T ) . halo 'C' in 
their labelling system. These halos were extracted from 
a cosmological simulation in a cube of comoving volume 
(lOOMpc) 3 and were chosen to have masses close to that of 
the Milky Way (~ 1O 12 M0) and avoid dense environments 
no neighbour excee ding half its mass within l/i _1 Mpc) 
Navarro et alj|2010h . 

As in Aquarius, we employ a 'zoom' resimulation tech- 
nique, with higher mass boundary particles used to model 
the large scale potential and lower mass particles in a 
~ 5/l _1 Mpc region surrounding the target halo. Extra power 
is added to the initial particle distribution on small scales 
in the high resolution region, as described by iFrenk et al.1 
(1996). We assume a ACDM cosmology, with parameters 
n m = 0.25, Ov = 0.75, tt b = 0.045, cr 8 = 0.9, n s = 1 and 
H = lOO/ikms^Mpc- 1 = 73kms- 1 Mpc _1 . 

The highest resolution realisation of halo C in Aquarius 
had a dark matter particle mass of 1.4 x 1O 4 M . However, 
the extra computational time associated with hydrodynamic 
simulations makes such a resolution impractical; our high- 
est resolution instead corresponds to a dark matter particle 
mass of ~ 2.6 x 10 5 Mq and an initial gas particle mass of 
5.8 x 10 4 Mq. In order to conduct convergence studies, we 
also simulated the halo at two lower resolutions, with par- 
ticle masses ~ 8 and ~ 64 times large r. We a dopt the same 
naming convention as ISpringel et al.l (|2008l ). labelling the 
three runs (in order of decreasing resolution) Aq-C-4, Aq- 
C-5 and Aq-C-6. Table [1] lists the numerical parameters of 
each simulation. 

Our simulation code is based on an early version of the 
PM-Tree-SPH code G ADGET-3. Baryonic pro cesses are mod- 
elled as described in lOkamoto et al.l (|2010l ) , with a num- 
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ber of modifications designed to improve the treatment of 
supernovae-driven winds. In the following subsections, we 
summarise some of the most important features of the code 
with emphasis on aspects that have the greatest impact on 
satellite formation. 



2.1 Radiative cooling and the equation of state 

Radiative processes in our mod el are implemented as de- 
scribed in W iersma et al.l l|2009al ) and include inverse Comp- 
ton scattering of CMB photons, thermal Bremsstrahlung, 
atomic line cooling and photoionisation heating from Hy- 
drogen and Helium. All gas in the simulation volume is 
ionised and heated by a spatiall y uniform, time evolvin g 
UV background, as calculated bv lHaardt fe Madaul l|200lf ). 
During the reionisation of H and He I (z = 9) and He II 
(z — 3.5), an extra two eV per atom of thermal energy is 
added to the gas, smoothed over Gaussian distributions with 
widths Az — 0.0001 and 0.5 respectively, in order to approx- 
imately account for non-equili brium and radiative-transfer 
effects (|Abel fe Haehneltlll999h . The contributions to heat- 
ing and cooling from eleven elements (H, He, C, N, O, Ne, 
Mg, Si, S, Ca and Fe) are inte rpolated from tables output by 
CLOUDY l|Ferland et al.| [T998). using elemental abundances 
smoothed over the SPH kernel, to approximate the mixing 
of metals in the interstellar medium (ISM). This avoids un- 
physical small scale fluctuations in the cooling time that 
arise if t he abundances associat ed with individual particles 
arc used l|Wiersma et al.ll2009bl l. 

Failure to resolve the Jeans mass (Mj) or Jeans 
Length (Aj) is known to lead to spurious fragmenta- 
tion through gravitat ional instability ([Bate fe Burkertll 19971 : 
iTruelove et al.l [l997h . Techniques employed to avoid this 
problem typically involve some form of energy injection 
to maintain an effective pressure th at guarantees that the 
available resolution is sufficie n t ("e.g.. iMachacek et al.l 2001; 
Robertson fe Kravtsovl |200S|; ISpringel fe Hernquistl bo03 ; 



Ceverino et al.ll2010l : ISchave et al.ll2oTol v" This pressure sup- 
port has been identified with, for exa mple, a hot phase main- 
taine d through energy input by SNe (|Springel fe Hernquistl 
2003) and tu rbulence induced by th e disk's self gravity 
and rotation (|Wada fe Normanl 120071 ). We include a min- 
imum pressure explicitly by adopting a polytropic equa- 
tion of state (EoS) with P m i n oc p 7cff for gas above the 
density threshold for star formation and adopt 7 e g = 1.4. 
For 7 e ff > 4/3, the Jeans mass increases with density 
(e.g.. ISchave fe Dalla Vecchial 2008). such that, if it is re- 
solved at the thres hold density, it is resolved everywhere. 
ISchave et af] (|2010l ) showed that, in models where star for- 
mation is strongly regulated by stellar feedback, as it is in 
ours, the choice of 7 e ff has very little impact on the global 
star formation rate. The star formation threshold density is 
Tin > 0.1cm -3 for the Aq-C-6 simulation and a factor of four 
and sixteen higher for the Aq-C-5 and Aq-C-4 simulations 
respectively. This scaling is chosen since, for irradiated pri- 
mordial gas with an isothermal density profile, halving the 
gravitational softening will increase the maximum density 
that is resolved by a factor of four. 



2.2 Multiphase ISM and Star Formation 

At sufficiently high pressure s, the ISM is known to exis t 
in distinct phases. We follow ISpringel fe Hernouis 3 l|2003h . 
modelling gas above the threshold density using hybrid SPH 
particles, which are assumed to consist of a series of cold 
clouds in pressure equilibrium with a surrounding hot phase. 
The total mass in the cold phase can increase through ther- 
mal instability and decrease through star formation and 
cloud evaporation by SNe, but the mass spectrum of clouds 
is kept fixed as 



$(m) = — — oc m 
dm 



(1) 



We adopt a = 1.7, guided by obse rvations that sug- 
gest a plausible r ange o f 1.5 — 1.9 ( Solomon fe Rivolol 



19891; iFukui et al l 1200 ll; iHever et~aH 1200 ll ). We follow 
Samland fe Gerhard! l|20ol ) in assuming clouds to be spher- 
ical, with size at a fixed mass d etermined solely by the am- 
bient pressure l|Elmegreenll 19891 ): 
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r(m) 



pc 



= 190P 



1/2 



(2) 



where P4 



- and k is the Boltzmann constant. The 
dependence of each cloud's dynamical time on the effective 
pressure follows directly from Eqn. [2] 
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and we assume that the star formation rate in each cloud is 
inversely proportional to its dynamical time (td vn )' 



tdy 



(4) 



where c» is the star formation efficiency, which is set to 
reproduce the no rmalisation of the Kennicutt-Schmidt law 
ijKennicut j I1998T ). The total star formation rate for each 
SPH particle is obtained by integrating Eqn.|4]over all clouds 
deemed capable of supporting star formation (which we as- 
sume to be in the mass range 10 4 — 1O 6 M0). 

SPH particles spawn new collisionless star particles in a 
stochastic fashion, with a probability that depends on their 
star formation rate and on the mass in the cold phase. Each 
star p article represent s a single stellar population, forming 
with a lChabrierl (|2003h initial mass function (IMF). Energy, 
mass and metals are returned to the ISM by AGB stars, type 
la and type II SNe on timescales appropriate for the age and 
metallicity of the stella r population, w i th yie lds a nd stel- 
lar lif etimes taken from IPortinari et al.l () 1998h and iMarigol 
(2001). Rather than performing this calculation at every dy- 
namical timestep, which is very expensive computationally, 
we use coarser steps, chosen such that the timescales asso- 
ciated with type II and type la SNe can be adequately sam- 
pled. Initially, each step has length isM o /50, where tsM is 
the lifetime of an 8M© star 0. When the age of the stellar 
population exceeds tgM Q , the timesteps lengthen to lOOMyr, 
which is short enough to model the release of mass and en- 
ergy from type la SNe and from AGB stars. 



t&Mp. ~ 40 Myr at solar metallicity 
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2.3 Supernovae Winds 

The perennial problem with the distribution of stellar feed- 
back energy in cosmological hydrodynamical simulations has 
been that the star-forming gas that receives the energy is 
dense enough to radiate it away before it can have any dy- 
namical effect ijKatz et al.lll996T l. This is likely to be a conse- 
quence of the inability to resolve the detailed structure of the 
ISM (e.g., Id alia Vecchia fc Schavell2008l ; ICeverino fc Klypinl 
2009). We employ a commonly used technique to circum- 
vent this limitation, which consists of imparting kinetic en- 



propelled is chosen at random to be parallel or anti-parallel 
to the vector (vq —v)x a gra v where vq is the velocity of the 



gas particle before it receives feedback energy, 



is the 



ergy to gas particles directly (e.g., Navarro fc White] 1 19931 : 



iMihos fc Hernauistlll994l : ISpringel fc Hernauistll2003h . The 
velocity we choose to give gas particles that receive SNe en- 
ergy is motivated by observations that suggest that large 
scale outflows have velociti es that scale w ith the circular ve- 
locity of their host galaxies (|Martinll2005l) . As a proxy for the 
host halo's circular velocity, which is computationally expen- 
sive to calculate for each particle on-the-fly, we use the local 
one-dimensional velocity dispersion, determined from neigh- 
bouring dark matter particles. This quantity is strongly cor- 
related with the maximum circular v elocity, ii mal , in a wa y 
that does not evolve with redshift (|Okamoto et al.| [2010). 
Our prescription results in a wind speed that increases as 
the halo grows and hence, from energy conservation, in a 
mass loading (wind mass per unit star formation rate) that 
is highest at early times. This scaling has been shown to 
give a much better match to the luminosity function of the 
Milky W ay satellites than mo dels that use a constant wind 
velocity jOkamoto et al.ll201oT ). 

One further addition to the model is needed 
to ensure that SNe drive n winds act as intended. 
iDalla Vecchia fc Schavel (2008) showed that standard kinetic 
feedback is more effective in low mass galaxies, where wind 
particles tend to drag neighbouring gas out with them. In 
high mass galaxies on the other hand, the pressure of the 
ISM can be sufficient to prevent much of the mass in the 
wind from escaping. Since we wish to be able to prescribe 
the mass loading and wind velocity directly, we choose to 
decouple wind particles from the hydrodynamic calculation 
for a short time in order to allow them to escape the high 
density star forming regions. When the density has fallen to 
Tin = 0.01cm -3 , the particles feel the usual hydrodynamic 
force again. If they do not reach sufficiently low densities 
after a time 10kpc/v w i n d, they are recoupled anyway. 

When a gas particle receives SNe energy from a neigh- 
bouring star particle, the wind speed (v m ) is obtained from 
the local velocity dispersion and then the particle is assigned 
a probability to be added to the wind: 



Pw 



AQ 



2 ITlsphVw 



(5) 



where AQ is the total feedback energy received by the gas 
particle and m sp h is the current mass of the SPH particle. 
Note that an SPH particle's mass may increase if it receives 
mass from SNe or AGB stars in neighbouring star particles, 
or decrease if it spawns a new star particle, which has a 
mass of half the original gas particle mass. If p m exceeds 
unity, that is, if there is energy available in excess of that 
needed to add the particle to the wind, then the extra energy 
is distributed to the gas particle's neighbours as an increase 
in internal energy. The direction in which wind particles are 



gravitational acceleration vector, pointing approximately to 
the local potential minimum (halo centre) and v is the bulk 
velocity of the halo, which we take to be the mean velocity 
of the gas particle's dark matter neighbours. The result of 
this treatment is a wind launched preferentially along an ob- 
ject's rotation axis (|Springel fc Hernquis t 2003). Our model 
for SN e winds differs from that described by Okamoto et alj 
|2010l) in two ways. Firstly, we allow all gas particles, not 
just those above the star formation density threshold, to be 
added to the wind if they receive feedback energy. The orig- 
inal prescription can result in a variable wind mass loading 
depending on how well the star forming region is resolved. 
Secondly, only type II SNe contribute to the winds, type la 
SNe energy is added to the gas as thermal energy. 



2.4 Satellite Identification 

Galaxi es are identified usin g a version of the SUBFIND algo - 
rithm l|Springel et al.ll200lf ) adapted bv lDolag etahl (|200ST) . 
which identifies self-bound structures and includes the inter- 
nal energy of gas when computing particle binding energies. 
From the ~ 5Mpc high resolution region, we select all galax- 
ies within 280kpc of the centre of the most massive (cen- 
tral) galaxy. This distance was chosen to match the limiting 
magnitude of the comple teness-corrected satell ite luminosity 
function constructed bv lKoposov et all (|2008l ). The largest 
satellite in our Aq-C-4 run is resolved with about 1.5 x 10 5 
particles in total, ~ 3 x 10 4 of which are star particles. In 
the following, we consider all galaxies with more than ten 
star particles, which, taking into account the typical mass 
fraction lost through stellar evolution for our choice of IMF, 
implies a stellar mass limit of ~ 1.2 x 10 5 Mq for Aq-C-4. 



3 THE EFFECT OF BARYONS ON 

SATELLITE DARK MATTER HALOS 

Using a dark matter only (DMO) counterpart of our 
Aq-C-4 run, simula ted as part of the Aquarius project 
i|Springel et al.ll2008h . we have examined the extent to which 
the dynamics of the baryons alter the structure of dark mat- 
ter (sub)halos of satellite galaxies over the course of their 
formation. The DMO run had identical initial conditions 
to our Aq-C-4, but for the absence of baryons and a cor- 
respondingly higher dark matter particle mass by a factor 

~ 1/(1 - Slb/tom). 

Naively, one might simply compare each subhalo with 
its DMO equivalent at z — 0, but this turns out to be prob- 
lematic. As has been noted in previous N-body simulations 
at different resolutions, small phase deviations in subhalo 
orbits get amplified over time, such that subhalos can be in 
quite different posit ions at z = (e.g.. iFrenk et al1ll999l ; 
ISpringel et ai] [2008). We see similar differences between Aq- 
C-4 and the DMO run. Subhalo orbits are also affected by 
other factors such as subhalo-subhalo scattering and varia- 
tions in the potential due to small differences in the growth 
history of the main halo. 

Since the strength of tidal s hocking is strongly 
dependent on pericentric distance (|Gnedin et alj 1 19991 : 
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iMaver et ai1l200ll ). small orbital deviations can cause large 
differences in subhalo structure, which are entirely unrelated 
to the presence or absence of baryons. This complication can 
be avoided, either by choosing subhalos with no close peri- 
centre, or by making the comparison at the epoch when the 
satellite is first accreted into the halo of the main galaxy, be- 
fore the orbits have had a chance to diverge. We choose the 
latter option, since the former restricts us to a very small 
number of cases, although we note that one massive halo 
in a low eccentricity (~ 0.2) orbit with a distant pericentre 
(~ 200kpc) shows comparatively small differences in its dark 
matter density profile at z = relative to the DMO case. In 
the few instances where the accretion times of the subhalo 
differ slightly between the hydrodynamical and DMO runs, 
we choose the earlier of the two epochs. 

In Fig. [T] we show spherically averaged profiles for the 
dark matter density (plotted as pr 2 to emphasise small dif- 
ferences) of our most massive satellites, for Aq-C-4 and the 
DMO run at the output time wh en each satellite fi rst joins 
the main friend-of-friends (FOF; lDavis et al.lll985T ) group 
The differences in the subhalo density profiles with and with- 
out baryons clearly exceed the uncertainties associated with 
finite sampling, indicated by the error bars. They are also 
greater than, for example, the differences expected between 
dark mat ter realisations of the same subhalo at different res- 
olutions |Springel et al.| [2008). Some subhalos (e.g. top left 
and centre right panels) appear to have been largely unaf- 
fected, whilst others (e.g. top and bottom right panels) show 
more substantial changes of up to thirty percent in some ra- 
dial bins. In Section [5J we describe an extreme example of 
a subhalo which suffered much more extensive damage as a 
result of baryonic processes. In general, however, there does 
not seem to be any consistent trend for baryons to increase 
or decrease the central density of the dark matter. 

In Fig. [21 we show spherically averaged dark matter 
velocity dispersion profiles for the same selection of satellites 
in Aq-C-4 and the DMO run. The largest differences are seen 
in those subhalos that show the most change in their density 
profiles in Fig. [1] Once again, the differences are typically 
less than ten percent in any given radial bin, but as much as 
thirty percent in some instances, with no apparent trend for 
baryons to raise or lower the dark matter velocity dispersion. 

From these results, we conclude that the baryons have 
had a relatively small impact on the dark matter phase-space 
structure of the subhalos, with the important caveat that it 
is unclear whether such effects are limited by the resolution 
of our simulations. Another important factor contributing 
to this conclusion is the strength of our feedback, since it 
dictates how easily baryons are able to condense in the cen- 
tre of low mass halos and affect the dynamics of the dark 
matter. Models with much weaker feedback might achieve 
more pronounced differences than we see here, but as we 
will show in Section 15.11 such models typically overpredict 
the luminosities of satellite galaxies. 

Although some satellites are accreted at fairly high red- 
shift (see the labels on each panel in Fig. [!}, it is unlikely 
that the baryons would have an increased effect in the re- 



2 In practice, we require the subhalo to be counted in the main 
FOF group for two consecutive snapshots to avoid instances where 
subhalos are spuriously joined to the main group for a short time. 



maining time to z = 0. In all but the largest satellites, gas 
is lost fairly rapidly following accretion, as we will show in 
Section [O] 



4 CONVERGENCE OF SATELLITE 
PROPERTIES 

In this section we investigate the convergence of various key 
properties of our simulated satellite galaxies. We make the 
comparison on an object-to-object basis, matching up satel- 
lites between runs. As explained in Section [31 this cannot be 
accomplished simply by choosing subhalos that are spatially 
closest at the final time. Instead, we trace particles back to 
the initial conditions and match them spatially there. In the 
following analysis we consider the most massive satellites in 
Aq-C-4 for which resolved counterparts exist in Aq-C-5 and 
often also in Aq-C-6. 

4.1 Stellar Mass 

We begin by considering the total stellar mass in each satel- 
lite at z = 0. This is a function of the rate at which gas can 
cool onto the galaxy and the efficiency of star formation, 
dictated by the gas physics and feedback. As well as check- 
ing that our results do not depend on resolution, we com- 
pare with an indepe ndent modelling technique, presented by 
ICooper et al.l (|2010t ). hereafter CIO. They used the Aquarius 
simulations to track the formation of dark matter substruc- 
tures in six different halos and a version of the semi-analytic 
galaxy formation code GALFORM to compute the baryonic 
properties of satellite galaxies0 At each output time, the 
stellar mass formed since the last simulation snapshot is as- 
signed to some fraction of the most tightly gravitationally 
bound dark matter particles in the subhalo, providing spa- 
tial and kinematic information for the stars. The 'tagged' 
fraction was chosen to match the distribution of sizes (half- 
light radii) for Local Group satellites and also to produce 
results robust to changes in resolution. 

Fig. demonstrates that the stellar mass in each sub- 
halo agrees relatively well between the three resolutions, 
particularly for the most massive examples, although the 
difference is as large as a factor of six in one case. Some of 
this scatter (between different resolutions and between the 
semi-analytic and hydrodynamical realisations) is likely re- 
lated to the deviations in subhalo orbits between simulation 
runs described in Section [31 Small differences in pericentric 
distance and eccentricity can strongly affect the tidal field 
and hence the extent to which stars can be stripped from 
subhalos. Examples of subhalos that have lost more than 50 
percent of their peak stellar mass (excluding the effects of 
stellar evolution) in one or more of the runs are indicated in 
Fig- El by circled points. 

The semi-analytic prescription typically predicts a lower 
stellar mass in each subhalo, although the correlation indi- 
cates that the ranking of subhalos by stellar mass is similar 

3 Their semi-ana lytic model is essentially that presented by 
iBower et al. | 1120061) . but with a lower circular velocity threshold 
(30 kms -1 ) to identify halos in which cooling is suppressed by 
reionisation. This value is motivated by recent hyd rodynamical 
simulations jHoeft et al.ll2006l ; [pkamoto et alj|2008h . 



© 201? RAS, MNRAS OOO.fTHTSI 



6 0. H. Parry et al. 



10" 



10" 



o 

Q_ 

i) 10" 



o 

"o 

* — i 

w 

5 

£ 10" 3 



10" 




10" 



Figure 1. The spherically averaged dark matter density profiles of the nine most massive satellite galaxies in our high resolution 
hydrodynamic (red) and dark matter only (black) runs. Error bars are obtained by bootstrap resampling. The comparisons are made at 
the redshift where the galaxy is first accreted as a satellite, which is shown as a label at the top of each panel. Black dashed vertical lines 
indicate the scale on which softened gravitational forces become fully Newtonian. Blue dashed vertical lines indicate the convergence 
radius of I Power et al.[ (2003). 



to that in our simulations. This offset between the two tech- 
niques is not obviously attributable to a single aspect of 
either model, but one mechanism that may be important 
is ram pressure stripping. In the semi-analytic model, any 
hot gas is instantaneously stripped from satellites upon in- 
fall. In combination with strong SNe feedback, this quenches 
star formation in satellites very rapidly following accretion. 
As we will show in the next subsection, ram pressure and 
feedback act to produce a similar effect in our simulations, 
but star formation is able to continue for an appreciable time 
after accretion and right up to z = for the most massive 
satellites. 



4.2 Mass evolution 

We now examine how satellites in the three runs acquire 
their dark and baryonic mass and form stars. Fig. [4] shows 
the gravitationally bound mass of dark matter, stars and gas 



for the nine most massive satellites as a function of redshift. 
For reference, we also include Tabled which lists the masses 
at z = 0. Reionisation (z — 9) is marked with a vertical 
dashed line and arrows indicate the accretion time, that is, 
the time when the galaxy first becomes a satellite in the 
high resolution run. Where the accretion times differ slightly 
between runs, the value for the high resolution case is shown. 

Apart from numerical convergence, which we will dis- 
cuss next, there are a number of interesting features to note 
in Fig. [ 31 many of which wer e also observed in the simula - 
tion s oflOkamoto et all |20ld ). IWadepuhl fc Springell (|2010l ) 
and ISawala et aL (2011). There are instances of satellites 
being periodically stripped of mass as they pass through 
pericentre, most obviously the satellite tracked in the bot- 
tom right panel, which loses dark matter and gas from its 
outer parts in two close approaches. None of these massive 
satellites appear to be on orbits with sufficiently high ec- 
centricity and/or a close pericentre to strip the more tightly 
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Figure 2. The spherically averaged ID velocity dispersion profiles of the most massive satellite galaxies in our high resolution hydro- 
dynamic (red) and dark matter only (black) runs. Error bars are obtained by bootstrap resampling. The comparisons are made at the 
redshift when the galaxy is first accreted as a satellite, which is indicated on the corresponding panels in Fig. [T] Black dashed vertical 
lines indic ate the scale on wh ich softened gravitational forces become fully Newtonian. Blue dashed vertical lines indicate the convergence 
radius of I Power et al. I j2003h . 



bound stellar component, although there are several exam- 
ples of satellites in the simulation that were heavily stripped 
or disrupted entirely and hence are not among the most mas- 
sive at 2 = 0. 

The loss of gas following accretion in most cases is fairly 
rapid and is brought about through a combination of ram 
pressure stripping and s tellar feedback. As noted by bot h 
lOkamoto fc Frenkl <|2009h and lWadepuhl fc Springel feoich . 
reionisation appears to have virtually no impact on satellites 
this large (the lowest mass example in Fig. [4] has final total 
mass of 5 x 10 8 Mq); the gas mass rises steadily through 
2 = 9, and star formation continues unabated. We note that 
the three largest satellites retain a substantial amount of gas 
and are still increasing their stellar mass at the present day, 
analogously to the ongoing star formation in the MW dlrrs. 

In terms of numerical convergence, in most respects, 
there is good agreement between the three resolutions. The 
exceptions to this are the rate at which gas is lost from 



some satellites following accretion and the resulting effect 
on the late-time star formation rates. In most cases, there 
is a clear tendency for more efficient ram pressure strip- 
ping with decreasing resolution. This affect appears to be 
related to the force resolution, which results in gas particles 
in lower resolution runs being less tightly bound and hence 
more susceptible to ram pressure stripping. The slightly dif- 
ferent timescales over which satellites are able to retain their 
gas and continue to form stars account, at least in part, for 
the often lower final stellar masses in Aq-C-5 and Aq-C-6 
noted in the previous subsection. 



5 OBSERVED PROPERTIES 

In the previous section, we demonstrated that our model 
produces satellite galaxies with properties that show rea- 
sonable convergence with resolution and stellar masses that 
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Figure 4. The mass evolution in dark matter (solid lines), stars (dashed lines) and gas (dotted lines) for the nine most massive satellite 
galaxies in our Aq-C-4 (red), Aq-C-5 (green) and Aq-C-6 (blue) simulations. Note that no resolved counterpart was identified in the 
Aq-C-6 simulation for the satellites tracked in the bottom row of panels. Stellar masses are scaled up by a factor of ten to reduce the 
range of the ordinate axis. The vertical dashed line indicates the redshift of reionisation and arrows the epoch at which the galaxy was 
first accreted as a satellite in Aq-C-4. The panels are ordered by total mass, from left to right along each row, beginning at the top. 
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Table 2. The mass in dark matter, gas and stars gravitationally 
bound to the nine most massive satellites at z = 

scale with subhalo mass in a fashion expected from an al- 
ternative modelling technique. We now proceed to exam- 



ine how well their observable properties match those of 
the Local Group satellites. Where photometric quantities 
are require d, we use the stellar population synthesis model 
PEGASS (|Fioc fc Rocca-Volmerangelll997l ). summing the 
luminosities of all star particles gravitationally bound to the 
subhalo at z = 0. 

5.1 Satellite Luminosity Function 

One of the most fundamental properties of any galaxy pop- 
ulation is its luminosity function. Encoded in its shape and 
normalisation are a range of physical processes that are key 
to understanding the formation and evolution of the popu- 
lation. 

In Fig. [5] we plot the luminosity function of sim- 
ulated satellites at each resolution. Also plotted is the 
(completeness-corrected) average luminosity function for 
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Figure 3. A comparison of the stellar mass that forms in each 
satellite at the three different resolutions of our hydrodynamical 
ru ns, as well as in the in dependent semi-analytic model presented 
by ICooper et al Note that not all satellites in the high 

resolution run have resolved counterparts at intermediate and low 
resolution. Points that lie on the abscissa correspond to subhalos 
with M sto n ar < 5 X 1O 4 M0. Data points are plotted for all satel- 
lites in the high resolution run with more than ten star particles. 
If the satellite also has ten star particles in the lower resolution 
run, it is plotted as a filled, rather than open symbol. Satellites 
that have lost more than 50 percent of their maximum stellar 
mass through tidal stripping are indicated by circled points. The 
dark and light gray shaded regions represent factors of two and 
five respectively away from the line of equality. 

MW and M31 satellites from lKoposov et al.l l|2008l ) that in- 
cludes the SDSS ultra-faint objects. We find good agreement 
between resolutions at all resolved luminosities, consistent 
with the convergence of stellar masses and star formation 
histories demonstrated in Section [4] The simulated popula- 
tions provide good matches to the Local G roup average at 
the fa int end, consistent with the findings of lOkamoto et alj 
(2010) in another of the Aquarius halos, but have no galaxy 
as bright as the Large Magellanic Cloud (LMC). The bright- 
est galaxy (in Aq-C-4) has a similar V-band magnitude to 
that of the Small Magellanic Cloud (SMC). 

Note that there is still a significant observational un- 
certainty in the total mass of the M W. Current estimates 
put it betw een 0.8 and 3 x 10 12 Mm dDehnen et ail 120061 ; 
iLi fc Whitell2008l ; IXue et al.ll2008l ; IWatkins et al.ll2010h . Our 
halo has a mass slightly closer to the lower end of this range 
Of M c rit,200 = 1.42 x 10 12 M Q . Clearly, if we are simulating 
a halo twice or half as massive as the MW, we should not 
expect to reproduce the luminosity function of its satellites 
exactly. A further consideration is whether the MW has a 
typical satellite population for its halo mass or total lumi- 
nosity. To this end, we have also plotted a black solid line 
and two shaded regions indicating the mean and spread of 
the luminosity function (see figure caption for details) for 
satellite systems in the SDSS around cen tral galaxies wit h 
r-band luminosity close to that of the MW (|Cuo et alj|201lf ). 
These data suggest that the Local Group is fairly typical, 
although satellites like the Magellanic Clouds are found in 
fewer than half of the systems in the sample. In another 



Figure 5. The rest frame, V-band luminosity function of satellite 
galaxies in our high (red), intermediate (green) and low (blue) res- 
olution simulations. Circles with error bars indicate the average 
MW+M 31 satellite lum i nosity function, corrected for complete- 
ness, by IKqdosov et all l|2008h . For consistency with their data, 
only simulated satellites within 280kpc of the central galaxy are 
included. The solid black line and shaded regions indicate the 
mean and spread in the satellite luminosity function s of Milky 
Way mass galaxies in an analysis of SDSS data by iGuo et al] 
(2011). The light grey region shows the RMS scatter about the 
mean in each bin. 

study, iLiu et all (|2010J ) found that > 80 percent of MW- 
like galaxies have no satellite as bright as the SMC within 
150 kpc. In this statistical context, perhaps the lack of very 
bright satellites in our simulations should not be a major 
cause for concern. 



5.2 Sizes 

An important and readily observable property of the most 
massive Local Group satellites is the distribution of their 
sizes, usually measured as the radius containing half the lu- 
minosity in projection. Unfortunately, a combination of the 
spatial resolution and the limitations of our subgrid model 
for star formation mean that we cannot hope to reproduce 
the observed sizes in our simulations. As explained in Sec- 
tion [2] a minimum pressure is maintained in star forming 
gas to ensure that the Jeans length on the equation of state, 
Aj,boS, is always resolved. Our subgrid model assumes that 
stars form on much smaller (unresolved) scales, inside molec- 
ular clouds, but the star particles that are created must 
nonetheless inherit the dynamical properties of the SPH par- 
ticle from which they formed. As a result, the minimum size 
of star-forming regions will be dictated by the warm/hot 
phase density and temperature, through \j,eoS, or by the 
gravitational softening, if this is larger. 

These limits are apparent in Fig. |JJJ where we plot 
the absolute V-band magnitude of simulated satellites as 
a function of their half-light radius. The observed half-light 
radii are de-projected by multiplying by a factor of 4/3, 
an approximation that is accurate to 2% for the exponen- 
tial, Gaussian, King, Plummer and Sersic profiles commonly 
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Figure 6. Absolute V-band magnitude as a function of de- 
projected half-light radius for satellite galaxies in our Aq-C-4 
(red filled circles) and Aq-C-5 (green filled circles) runs. The 
Magellanic Clouds are represented by purple crosses. Their pro- 
.jecte d half-light radii were deriv ed from angular values taken 
from [B othun fc Thompson! jl98S|) and distance e stimates from 
iHilditch et al] j2005h and IPietrzvriski et al] j200Sft . 3D half-light 
radii wer e calculated by m ultiplying these values by 4/3, as sug- 
gested bv lWolf etldl ||2010| '). All other MW satellites brighter than 
—7.5 (except Saggitarius) are shown as blue triangl es with lu- 
minos ity and de-projected half-light radii taken from IWolf et al] 
The solid blue line is a least-squares fit to all the observa- 
tional data points. Arrows indicate the scales on which softened 
gravitational forces become fully Newtonian for each resolution. 

used to fit the MW satellites l|Wolf et alj|201Ch . The stars 
in both the Aq-C-4 and Aq-C-5 runs typically have much 
less concentrated distributions than the observed satellites. 
An exception to this is the third brightest satellite in Aq-C- 
4, which has a V-band magnitude of —12.2 and a half-light 
radius of ~ 480pc. It has a very high mass fraction in stars 
and an unusual history, forming in a series of violent major 
mergers at z ~ 4 before being subjected to strong tidal dis- 
ruption between z = 2 and 2 = 0. We discuss this satellite 
in detail in Section [6] 

In the highest resolution run, Aq-C-4, we expect the 
gravitational softening to be the main factor limiting the 
minimum sizes of star forming regions, since it is always 
larger than Xj.eoS- In Aq-C-5 and Aq-C-6 (not shown here), 
which have lower threshold densities for star formation by 
factors of four and sixteen respectively, Xj,eoS at the thresh- 
old is comparable to the softening, so should also be impor- 
tant in setting the sizes of the stellar component. For both 
Aq-C-4 and Aq-C-5, the half-light radius of the most mas- 
sive galaxy should not be limited by either effect and is 
consistent with the observations, given the large scatter. 

5.3 Dynamical Masses 

The stellar kinematical properties of Local Group dwarf 
galaxies provide an important test of the ACDM cosmol- 
ogy. Subhalos that form in N-body simulations of MW-mass 
systems appear to have potentials compatible with the stel- 
lar k i nematics of the brig htest MW satellites l)Stoehr et al.l 
20021; IStrigari et alJl2010h . Nonetheless, the analytic calcu- 
lations required to reach such conclusions necessarily in- 



Figure 7. The total mass contained within the projected half- 
light radius (Rhi) of simulated satellites, measured in the high 
resolution, dark matter only, Aquarius-C-2 simulation (open and 
filled black circles). The radii are inferred from the observed 
luminosity-size relation (see text for details). The filled circles 
indicate those satellites that are used in the statistical compari- 
son with the observed data, illustrated in Fig. [8] The error bars on 
the simulated data points show the masses obtained by assuming 
half-light radii lc above and below the mean fitted values, where 
a is determined by the scatter about the fit to the observed data . 
Red crosses with error bars are estimates from[yVolf et al. d2010T ) 
for a selection of MW satellites. 

elude simplifying assumptions. Hydrodynamic simulations 
attempting to model star formation self consistently in a cos- 
mological setting are inevitably some way behind the best 
N-body simulations in terms of resolution and must also 
model uncertain baryonic physics on sub-kiloparsec scales. 
As such, our simulations are not suitable for studying the de- 
tailed kinematics of the stars directly; instead, we resort to a 
somewhat cruder comparison and ask whether our satellites 
form in realistic potential wells, by comparing simulated and 
observationally determined masses. 

Historically there has been significant uncertainty asso- 
ciated with determining satellite masses from observations. 
Typically, estimates are derived from the line-of-sight stellar 
velocity dispersion with three key assumptions: i) the system 
is spherically symmetric, ii) stellar orbits are isotropic and 
iii) the system is in equilibrium. Two recent studies have at- 
tempted a more general approach, with the aim of reducing 
the systematic uncertainties. Using an app r oach based on 
the spherical Jeans equation, IWalker et al] l|2009l ) showed 
that for the brightest MW dSphs, the mass within the pro- 
jected half-light radius is robust to changes in the anisotropy 
and underlying density profile. T his relation was explained 
analytically by IWolf et al] (|2010T l who demonstrated that, 
if the stellar velocity dispersion profile rem ains relatively 
flat in the centre, as observations suggest (e.g. JWalker et al] 
120071 ). then the uncertainty introduced by assuming a par- 
ticular anisotropy is minimised at the (3D) radius where 
the logarithmic slope of the stellar number density profile, 
— dlnn./dlnr- = 3. They also showed that, for a range of re- 
alistic light profiles that have been used to model the MW 
dSphs, this minimum lies close to the (de-projected) half- 
light radius. It is this radius, therefore, at which we choose 
to compare the enclosed masses of satellites in the simula- 
tions and observations. 
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In the previous subsection we described how aspects of 
our simulations, particularly the limitations of the subgrid 
treatment of the ISM and the gravitational softening scale, 
can set an artificial lower limit to the sizes of the stellar 
components of the satellites. However, we also demonstrated 
that the luminosity function of the simulated satellites is 
close to that observed, the stellar mass in each satellite is 
relatively well converged and stellar mass is found to scale 
with subhalo mass similarly using an alternative modelling 
technique. 

With these checks in mind, we proceed with the as- 
sumption that the cooling, star formation and feedback pre- 
scriptions in our model result in a realistic stellar mass in 
each satellite, but that stars form in a configuration that is 
too diffuse. We then ask what the projected half-light radius 
of each simulated satellite should be at a fixed luminosity, 
based on the observed sizes of the brightest MW satellites. 
We take a simple least-squares fit to the data points in Fig. [5] 
(minimising the sum of the squared differences in the mag- 
nitude coordinate) and compute the scatter in the (log) ra- 
dius coordinate about this line. For each simulated satellite, 
we assume a Gaussian distribution of possible sizes, with a 
mean equal to the fit evaluated at the satellite's luminosity 
and dispersion defined by the observed scatter. 

Following this procedure we find that, for magnitudes 
Mv> —12, our fit to the observed data implies sizes be- 
low the scale at which softened gravitational forces become 
non-Newtonian in the Aq-C-4 simulation, which leads to an 
underestimation of the enclosed mass. We choose instead 
to measure the mass of each satellite in a much higher 
resolution dark matter only realisation of the simulation, 
Aquarius-C-2, which has a smaller softening scale by a fac- 
tor of ~ 4, such that the fitted half-light radii of satellites 
down to My~ — 7.5 are larger than the force resolution. 
The central masses measured in the higher resolution sim- 
ulation are typically forty to eighty percent higher for the 
ten brightest satellites, but the difference can be a factor 
of three for satellites with Mv~ —8. Given the results pre- 
sented in Section[3] we do not expect the omission of baryons 
from Aquarius-C-2 to have had a large impact on the central 
densities and hence the measured masses of these satellites. 

Fig. [7] shows the mass enclosed within the mean half- 
light radius chosen for each satellite, with error bars indi- 
cating the masses corresponding to ±ler sizes. Although the 
range of plausible values is large, the brightest simulated 
satellites have mean masses three to five times higher than 
the MW satellites of the same luminosity. While there is less 
of a discrepancy at fainter magnitudes, our model seems to 
show a more gradual increase in luminosity with mass than 
is suggested by the observational data. We note that, start- 
ing from identical initial conditio ns to our Aq-C-4 simula- 
tion, [Wadcpuh l fc Springell l|201(]| ) found that the mass-to- 
light ratios of their satellites were typically higher than those 
quoted observationally by a very similar factor and were also 
more discrepant in the most massive satellites (see their Fig. 
15). In a hydrodynamic simulation of th e Local Group , again 
with resolution similar to our Aq-C-4, iKnebe et all (|2010h 
found a similar result for satellites bound to their MW and 
M31 analogues, with mass-to-light ratios a factor ~ 7 too 
high. We note, however, that they measured half-light radii 
for their satellites using the star particles forming in their 
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Figure 8. The cumulative fraction of the 4th — Ylth brightest 
satellites as a function of the mass contained within the half- 
light radius in three model distributions (solid b lue, green and red 
lines) and those derived for MW satellites by IWolf et~aH j201Ch 
(dashed line) . Labels in the top left indicate the number of times 
we repeat the process of drawing masses for our sample of nine 
satellites in order to define the model distribution, along with 
the probability that each model is consistent with the observed 
masses. 

simulation, which, as we have shown, can be too large when 
the scales associated with star formation are not resolved. 

To quantify the discrepancy in Fig. [7] statistically, we 
construct multiple realisations of the half-light masses of the 
simulated satellites by drawing sizes from the distributions 
described above and computing the mass enclosed in the 
corresponding high resolution Aquarius-C-2 satellites. We 
then combine the samples to define a model distribution for 
the cumulative fraction of satellites with mass larger than a 
given value and calculate the probability that the observed 
masses could have been drawn from it, using a one-tailed 
Kolmogorov-Smirnov (KS) test. 

In order to make a like- f or-lik e comparison with the 
masses quoted by IWolf et "all (|2010T ) , we consider only the 
4 th - 12 th brightest simulated satellites in our Aq-C-4 sim- 
ulation that are within 280kpc, corresponding to all of the 
classical MW satellites except the Magellanic Clouds and 
Sagittarius, and including Canes Venatici (which is approx- 
imately the same luminosity as Draco). The sample is indi- 
cated in Fig. [7]by filled circles. Note that the third brightest 
satellite in Aq-C-4 is in the midst of tidal disruption (see Sec- 
tion [6]), a process that, as a result of small differences in the 
orbits of the subhalos, is already complet e in Aquari u s-C-2 , 
hence no counterpart is found. Although IWolf et al.l (|2010l ) 
derived masses for fainter satellites, these correspond to sim- 
ulated galaxies with fewer than ten star particles, whose lu- 
minosities are uncertain in our simulations and which we 
therefore choose to exclude. In Fig. [8] we plot cumulative 
model distributions, drawing masses for the sample of nine 
satellites multiple times to define each distribution, as indi- 
cated by the labels in the top left of the plot. The probabil- 
ities that the observationally derived masses are consistent 
with those distributions are found to be around six percent. 

As both Fig.[7]and Fig.[8]demonstrate, the masses of the 
brightest simulated satellites are too high compared to those 
derived for the MW satellites. This is another manifestation 
of the problem recently highlighted bv lBovlan-Kolchin et al.l 
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(2011) who compared the measured masses within the half- 
light radii of the same satellites considered here with results 
from high-resolution simulations of cold dark matter halos, 
includi ng the Aquarius suite . Assu ming an NFW density 
profile (|Navarro et al.lll996bl . Il997l ), they showed that the 
most massive subhalos in the simulations are too concen- 
trated to be able to host the brightest observed satellites. 

The mismatch seen in Fig. [7] and in the results of 
iBovlan-Kolchin et al.l (|201ll ) could, in principle, be due to an 
underestimate of the central masses of t he observed satel - 
lites. However, for the errors quoted by IWolf et all (I2010T ) 
to be substantially underestimated would require rather ex- 
treme variations in the anisotropy profile, which would be 
poorly fit by their fairly general parameterised form. This 
seems unlikely to be the sole source of the disagreement be- 
tween model and data. 

While the discrepancy could be simply due to statis- 
tics, it might also reflect a serious shortcoming either of the 
standard CDM cosmogony or of current models of galaxy 
formation, such as those assumed in our simulations. A pos- 
sible explanation of the discrepancy between the mass-to- 
light ratios measured for the real and simulated satellites 
is that the central dark matter densities predicted in the 
CDM model are reduced by baryon ic physics. One mecha - 
nism for achieving this, proposed bv lNavarro et al.l (|l996al ). 
is the condensation of a dense baryonic component followed 
by the rapid expulsion of gas by stellar feedback. The dark 
matter adjusts to this change in the potential by develop- 
ing a central "core", shifting the rotation curve maximum 
to a larger radius and reducing the mass-to-light ratio in 
the central parts. This process does indeed appear to play 
an important role in the evolution of one satellite in Aq-C- 
4 (see Section [6]), which forms in the subhalo that has the 
largest mass prior to accretion. If this process is common, it 
is possible that it is not seen here in less massive subhalos 
due to lack of resolution. 

A more radical explanation of the discrepancy is that 
the dark matter consists of warm, rather than cold, par- 
ticles. In this case, subhalos of a given mass form later 
and have lower con c entrations than in the CD M model (see 
iNavarro et~alll 19971 ; iHogan fc Dalcantonllioool ') . lLovell et ail 
1 20 111 ) have recently shown explicitly that the masses and 
concentrations of subhalos in a warm dark matter model 
agree well with the data. 



6 A STAR-DOMINATED SATELLITE 

The formation history of one of the satellite galaxies in our 
high resolution hydrodynamical simulation is particularly 
interesting. By z = 0, we find that it has become dominated 
by its stellar component, with a mass-to-light ratio of ~ 2.4 
and its dark matter has become much less concentrated than 
otherwise similar subhalos. It appears as an outlier in Fig- EI 
as it has a very small half-light radius for its luminosity. In 
this section, we briefly describe its formation history and 
explain why it develops into such an unusual object. 

At z = the satellite is, in fact, in the process of being 
tidally disrupted and has a substantial stellar stream asso- 
ciated with it. Fig. illustrates the structure of the stream 
in two orthogonal projections centred on the main galaxy. 
The dense stellar nucleus of the satellite that remains iden- 



tifiable as a bound structure is also visible. We track all star 
particles associated with the satellite at the epoch when it 
is accreted and plot their projected mass density at z = 0. 
The stripped stars account for the majority of the stellar 
halo by mass. 

The stream is a result of a fairly eccentric orbit with 
several close pericentres, illustrated in the bottom panel of 
Fig. 1101 which shows the distance of the satellite from the 
centre of the main galaxy as a function of redshift. The 
dashed line indicates the virial radius of the main halo. The 
accretion time is the point where the two lines intersect. 
The top panel tracks the mass in gas, dark matter and stars 
bound to the satellite over the same redshift interval, as 
well as the total mass fraction in stars. At accretion, it is 
the brightest satellite of the central galaxy, but only the 
third brightest at z = 0, as a result of the reduction in 
stellar mass through tidal stripping. The stellar fraction at 
accretion (~ 0.02) is fairly typical of the surviving satellites. 
Note that it is very common for the stellar fraction of a 
satellite to increase with time after it is accreted, since the 
outer parts of the dark matter halo are less tightly bound 
than the stars and hence more susceptible to t idal stripping 
fe.g. lPenarrubia et al.ll200Sl ; ISawala et al.ll201ll ). The middle 
panel of Fig. [TO] shows the evolution of the central (r < lkpc) 
density of the subhalo in gas and dark matter, both of which 
drop sharply when the satellite is close to pericentre. A de- 
cline is also evident after z ~ 3, well before the satellite is 
accreted, the origin of which we discuss in more detail below. 

During the first few orbits, the stellar component re- 
mains unaffected while the dark matter lying beyond the 
radial extent of the stellar component is stripped. In fact, 
some of the dark matter particles with pericentres within 
the stellar component are also stripped, as a result of the 
dark matter having a higher radial velocity dispersion than 
the stars. The final masses of the stellar and dark matter 
components are factors of ~ 50 and 2 x 10 4 lower than their 
peak values respectively. 

The extent to which the two components are stripped is 
strongly affected by their radial density profiles, which are 
shown in Fig. [TTJ at the time of accretion. The overplotted 
regions indicate the range of densities (±lcr) in each bin 
for the nine most massive surviving satellites. Clearly, the 
stellar component of this galaxy is unsually concentrated 
relative to those other galaxies, whilst the dark matter and 
gas have shallower than average central density profiles. It 
is unclear how much effect the gravitational softening has in 
this respect, since forces begin to become sub-Newtonian on 
scales less than < 720pc, but we note that the profiles also 
differ outside that radius. The highly 'cusped' stellar profile 
allows the central stellar nucleus to resist the strong tidal 
forces that unbind the majority of the dark matter. It also 
accounts for the unusually small half-light radius shown in 
Fig.H 

The origin of these density profiles is related to the 
satellite's violent formation history. In a series of major 
mergers at z ~ 3, gas is funnelled to the centre of the main 
progenitor, initiating an intense burst of star formation that 
gives rise to a highly concentrated stellar distribution. The 
subsequent burst of feedback energy rapidly removes a large 
fraction of the gas and leads to a fall in the mean binding 
energy of the central dark matter. This episode is clearly vis- 
ible in the bottom panel of Fig. 1101 which shows the central 
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Figure 9. Two orthogonal projections of the stellar stream associated with the satellite. The centre of the main halo is at the origin. The 
stream is defined by selecting all stars associated with the satellite at the time of accretion and locating them at z = 0. The surviving 
satellite is clearly seen as the bright concentrated object lying along the stream. 



gas density drop as it is expelled by feedback and turned 
into stars, followed by a decline in the dark matter density 
in response to the change in the potential. This sequence 
of events is e ff ectivel y the process originally proposed by 
iNavarro et all dl996ah . The reduced binding energy of its 
central dark matter, along with an extreme orbit, combine 
to produce the unusual properties of this satellite at z = 0. 



7 CONCLUSIONS 

We have investigated the formation and evolution of a Milky 
Way-like satellite system in an SPH simulation over three 
levels of resolution, in which particle masses vary by a fac- 
tor of 64. The properties of our simulated satellites show 
relatively good numerical convergence, with the final stel- 
lar masses typically agreeing to within a factor of two, and 
always to within a factor of six. We also compared to an in- 
dependent estimate of the stellar mass expect ed to form in 
each s ubhalo, using the semi-analytic model of lCooper et alj 
(2010). The two theoretical techniques produce a similar 
ranking of the subhalos by stellar mass, although our simu- 
lations typically form a higher mass of stars by a factor of 
between two and six. This discrepancy may be partly ex- 
plained by the assumption in the semi-analytic model that 
gas is stripped instantaneously when a galaxy becomes a 
satellite. The mass evolution in gas, stars and dark matter 
of each satellite agrees well between resolutions, except that 
gas is stripped more rapidly at lower resolution following 
accretion onto the main halo. Poorer force resolution causes 
gas particles to be more loosely bound to the subhalo and 
hence more susceptible to ram pressure stripping. This phe- 
nomenon may account for many of the differences in the 
final stellar masses between resolutions. 

By comparing the dark matter halos of our satellite 
galaxies to those that form in a dissipationless version of 
the same simulation, we were able to quantify the expected 



impact of baryons on the phase-space structure of the dark 
matter. Due to small deviations in satellite orbits between 
different realisations of the same halo, it is necessary to make 
this comparison when the satellite first falls in rather than 
at z = 0. Although in some radial bins, in a few subha- 
los, the density and velocity dispersion profiles are found 
to change by ~ 30 percent, the differences were typically 
less than 10 percent. With the caveat that the resolution 
of our simulation may limit the magnitude of such effects, 
we conclude that baryons have a relatively small impact on 
the structure of the dark matter halos of satellite galaxies 
around MW-like hosts. 

Our model provides a reasonable match to the faint end 
of the Local Group satellite luminosity function averaged 
between the Milky Way and M31, although there is a slight 
deficit at th e bright end, with no LMC ana logue. However, 
SDSS data (|Liu et alJl201Ct iGuo et alj|201lh suggest that it 
is quite common for galaxies with luminosities like the MW 
to have no satellites as bright as the LMC and SMC. 

Due to the limitations of the spatial resolution and the 
implementation of baryonic physics in our simulations, par- 
ticularly the modelling of the multiphase ISM, stars do not 
form in sufficiently concentrated distributions to match the 
half-light radii of Local Group satellites. However, the rea- 
sonable agreement between the stellar masses in simulations 
with different resolution and different modelling techniques, 
combined with the match to the observed satellite luminos- 
ity function, suggest that the baryonic mass that is able to 
cool in each (sub)halo and form stars is realistic. 

In order to test whether satellites of a given luminosity 
form in halos with masses consistent with those observed, we 
compare their 'half-light masses' with val ues derived f or a se - 
lection of the brightest MW satellites by I Wolf et ail l|2010h . 
For this comparison, each satellite is assigned a distribu- 
tion of half-light radii from the best fit to the observed 
luminosity-size relation and its variance. In the hydrody- 
namical simulation, the gravitational softening is compara- 
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Figure 10. Top panel: the mass in gas (blue), dark matter (black) 
and stars (red) gravitationally bound to the satellite's main pro- 
genitor, and the stellar fraction (magenta, dashed line, measured 
on the right vertical axis) as a function of redshift. Some of the 
small variations in the masses (and hence also the stellar fraction) 
at z < 1 are due to the difficulty in identifying the subhalo's par- 
ticles against the high background density at the centre of the 
main halo. Centre panel: the density of dark matter (black) and 
gas (blue) within the central lkpc. Bottom panel: the distance to 
the satellite from the centre of the main halo. The dashed line 
indicates the virial radius of the main halo. 



ble to these fitted half-light radii, so we instead measure the 
masses in a much higher resolution, dark matter only re- 
alisation of the simulation. We hence explicitly ignore any 
effects baryons may have had on the central density profiles, 
which, in any case, the results in Section[3]suggest are small. 
The large scatter about the observed relation translates into 
a broad range of possible masses for each simulated satellite, 
but nonetheless, the mean masses for the brightest examples 
(Mv < — 11) are about three to five times higher than their 
observed counterparts. The observed mass-luminosity rela- 
tion seems to be somewhat steeper than that produced by 
our model, although due to the small sample sizes, the slope 
of the relation and the scatter about it are relatively poorly 
defined in both cases. A KS test, taking into account the un- 
certainties in the half-light radii assigned to the simulated 




Figure 11. The dark matter (black), stellar (red) and gas (blue) 
density profiles of the galaxy at z ~ 2, when it is first accreted as 
a satellite. Regions of the same colour indicate the spread of val- 
ues (±1<t) in each radial bin for the nine most massive surviving 
satellites at z = 0. 

satellites, returns a six percent probability that the observed 
masses could have been drawn from the distribution defined 
by the simulation data. 

Although the apparent disagreement between the simu- 
lations and the data could be simply due to statistics, there 
are also a number of plausible physical explanations. It could 
be that baryonic processes significantly reduce the central 
dark matter densities of satellite galaxies. Possible mecha- 
nisms to achieve this include, for instance, a sudden change 
in the local potential, induced by the rapi d expulsion of 
baryon ic mass through stellar feedback (e.g. iNavarro et al.l 
Il996al ) o r heating due to bulk m otions of dense clumps of 
gas (e.g. [Mashchcnko et al. 20061). The results in Section [3] 
imply that, if such processes are important, they are either 
not resolved in the Aq-C-4 simulation, or are not properly 
captured by our feedback prescription, except in one case, 
which happens to be the most massive subhalo at accretion. 
Less concentrated dark matter profiles would also result if 
the dark matter consis ts of warm, rather than cold particles 
(e.g. lLovell et alj|201lt) . 

The broad range of possible explanations of the discrep- 
ancy highlighted by our results illustrates how uncertain our 
understanding of galaxy formation still is on the scale of 
dwarf galaxies. Determining which, if any, is correct will be 
of critical importance in assessing the viability of the CDM 
cosmology and the success of galaxy formation models. 
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